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We analyze the electronic properties of a simple stacking defect in Bernal graphite. We show that 
a bound state forms, which disperses as \k — K\^ in the vicinity of either of the two inequivalent zone 
corners K. In the presence of a strong c-axis magnetic field, this bound state develops a Landau 
level structure which for low energies behaves as oc \nB\^^'^. We show that buried stacking 
faults have observable consequences for surface spectroscopy, and we discuss the implications for 
the three-dimensional quantum Hall effect (3DQHE). We also analyze the Landau level structure 
and chiral surface states of rhombohedral graphite, and show that, when doped, it should exhibit 
multiple 3DQHE plateaus at modest fields. 

PACS numbers: 81.05.Uw, 61.72.Nn , 73.43.-f 



I. INTRODUCTION 

An explosion of research activity associated with the 
novel two-dimensional material graphene has prompted 
a reexamination of its bulk parent, graphite. Much, of 
course, is known about graphite [Ij. Bernal graphite is a 
hexgonal crystal consisting of graphene sheets stacked in 
an ABAB configuration. The sp^-hybridized a electrons 
form double bonds between the carbon atoms, while the 
remaining tt electrons, in the orbital, are itinerant. 
The electronic structure parameters for graphite were 
first derived by Wallace, and by Slonczewski, Weiss, and 
McClure (SWMC) [Sj]. Within each plane, the tt electrons 
move on a honeycomb lattice with a nearest neighbor 
hopping integral 7q ~ 3.2 eV. Of the four atoms per unit 
cell, two are arranged in vertical chains, with a vertical 
nearest neighbor hopping of 7^ ^ —390 meV. Additional 
further neighbor hoppings are also present. For example, 
the TT electrons on the non-chain sites undergo two-layer 
vertical hopping through open hexagons in the neighbor- 
ing layers, with amplitude « — lOmeV. This results 
in a very narrow band of width 40meV along the K-H 
spine of the Brillouin zone, with electron pockets at K 
and hole pockets at H [4]. 

Recently, striking experimental observations of what 
may be bulk three-dimension al q uantum Hall plateaus in 
graphite has been reported [l^]- Any two-dimensional 
(2D) system, such as graphene, which exhibits the quan- 
tum Hall effect (QHE) should exhibit a 3DQHE if the 
interplane coupling is sufficiently weak. The reason for 
this is that the cyclotron gaps between Landau levels 
narrow continuously as one adiabatically switches on the 
c-axis couplings, and cannot collapse immediately. For a 
3D electron system in a periodic potential and subject to 
a magnetic field, a generalization of the TKNN result 0] 
by Halperin [6] shows that the conductivity tensor must 
be of the form 

whenever the Fermi level lies within a bulk gap, where 



eijk is the fully antisymmetric tensor and G is a reciprocal 
lattice vector of the potential (which may be G = 0). 
The Hall current is then carried by a sheath of chiral 
surface states. Eventually, however, the c-axis hopping 
will become large enough that the Landau gaps collapse. 
Equivalently, for a given value of the c-axis hopping 7^, 
the magnetic field B must exceed a critical strength 
in order that the Landau level spacing overwhelms the 
c-axis bandwidth and opens up a bulk gap. 

Typically, the field scale B^ is extremely large, and 
much beyond the scale of current experimentally avail- 
able fields. For a system with ballistic dispersion de- 
scribed by an effective mass m*, the orbital part of the 
spectrum {i.e. neglecting Zeeman coupling) yields a dis- 
persion = (n+^) huj^, where n is a nonnegative integer 
and where uo^ = eB/m*c is the cyclotron frequency. The 
cyclotron energy may be written as hw^ = W^^ • {B/B^) 

and the field scale as B^ = hc/ef2, where i? the unit 
cell area. W^^ is on the order of the bandwidth in zero 
field, which is typically several electron volts. Since 
hc/e = 4.13 X 10^ TA , B^ is typically enormous, on 
the order of tens of thousands of Tesla. Thus, if the 
c-axis bandwidth is W^, the critical field is given by 
B^ = (W^/W||) • 5^, and even highly anisotropic ma- 
terials with < 10~^ W^ii will have critical fields in the 
range of hundreds of Tesla. 

As shown by Bernevig et al. [5I] , similar considerations 
would apply for graphene sheets in AAA (simple hexag- 
onal) stacking. The Landau level dispersion is then 

E^{B, k) = 271 cos{k,c) + sgn(n) j^^nB/Bo , (2) 

with ^0 = B^/27rV3 = 7275 T, where 70 ^ 3.16 eV is 
the in-plane hopping and 7^ ~ 0.39 eV is the hopping 
integral between layers [3]. The gap between Landau 
levels n and n + 1 collapses at a critical field 

B..= M ^ (3) 
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For n = one finds B^q ^ 1800 T. However, due to the 
Bernal stacking, one finds [5|] that the principal cyclotron 
gap surrounding the n = Landau levels opens above 
Be = l^T (electrons; n = +1) or above Bc = 7T (holes; 
n = —1). When the Fermi level lies within either of these 
gaps, the Hall conductance is quantized at cr^y = 2e^ /hd^ 
where d = 3.37 A is the inter-plane separation. 

The analysis of ref. [5] shows that the second cyclotron 
gap should not open below fields on the order oi B^2 ^ 
1000 T. This suggests that the multiple QHE plateaus 
observed by Kempa et al. [l^ are of a different origin, 
and are not describable by a model of crystalline Bernal 
graphite alone. 

In this paper, we consider two variations which lead to 
a different plateau structure to that of crystalline Bernal 
graphite. The first is rhombohedral graphite, which is 
stacked in ABCABC fashion. For this structure, we find 



B. 



(Vn + 1 + ^/n^ 



5. 



c,0 



with o ^ 0.123 T. When lies in the Fermi level 
between the n and n + 1 Landau levels, the Hall con- 
ductivity is given by cr^y = (4n + 2)e^/hd. Ah initio 
calculations show that the total energy of rhombohedral 
graphite to be approximately 0.11 meV per atom larger 
than the Bernal hexagonal phase Q. With such a small 
energy difference, even highly oriented pyrolytic graphite 
(HOPG) is believed to contain several percent rhombo- 
hedral inclusions. Powdered graphite samples with up to 
~ 40% of the rhombohedral phase are obtainable [9]. 

The second possibility we examine is that of a sim- 
ple stacking fault in Bernal graphite, of the form 
ABABCBCB, This fault interpolates between two de- 
generate vacua - the ABAB and CBCB Bernal phases. 
We analyze the c-axis transport through such a defect, 
within a simple model of nearest neighbor hopping, and 
compute the S'-matrix as a function of in-plane wavevec- 
tor. As expected, the transmission is sharply attenu- 
ated in the vicinity of the Dirac points. We also find 
a novel bound state associated with the stacking defect, 
with two-dimensional dispersion E(k) ex \k — K\^ near 
the Dirac points. In the presence of a c-axis magnetic 
field, this leads to a bound state Landau level energy 
£^(n, B) (X |n5p/^. In the appendix, we undertake a cal- 
culation of the bound state spectrum in zero field for the 
full SWMC model [3] , which includes seven tight binding 
parameters. 

We conclude with a discussion of surface spectroscopy 
of buried stacking faults, and with remarks about the 
relevance of our results to future experiments. 



II. RHOMBOHEDRAL GRAPHITE 

In rhombohedral graphite (KG) there are two sublat- 
tices, in contrast to four in the case of Bernal hexagonal 




FIG. 1: Crystal structure of rhombohedral graphite. 



graphite (BHG). The primitive direct lattice vectors are 



a. 



a 



a X 



2 — 2<^^^+2^^ 

a^ = \ ax ^ ay ^ dz . 

The basis vector S = — + 0^2) separates the A and 
B sublattices. Note that = d z — 6. The lattice pa- 
rameters are a = 2.46 A and d = 3.37 A. 

Our treatment starts with a simplified version of the 
work of McClure [13] • We consider several types of hop- 
ping processes: 



(i) in-plane A — B hopping: 



AB 



-7o 



t{S) +t{a^+S) +t{a2 + S) 



(4) 



where, t{d) is a translation operator through a vec- 
tor d. 

(ii) neighboring plane diagonal A — B hopping: 

= 73 - as + ^) + ^(«2 - ^3 + ^) (5) 
+ t{ai + a2 — ttg + 5) 

(iii) nearest neighbor and next nearest neighbor plane 
vertical A — B hopping: 



(iv) neighboring plane diagonal A — A hopping: 

t{a^) -\- t{a^ — a-^) -\- t{a^ — + H.c. (7) 



ztAA _ 
^4 — 73 



(v) neighboring plane diagonal B — B hopping: 



H^^ t{a^) ^t{a^- a-^) ^t{a^- a^) + H.c. (8) 
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The full Hamiltonian is then given by 

Hi 



AB ^ jjAB 



H 



(9) 

where H^"^ = [H^^)^ for n = 1, 2, 3. From Wallace and 
SWMC H, we take 



7o 3160 meV 
72 = lOmeV 



7i = 390 meV 
73 250 meV . 



(In the language of McClure [10], 72 = 72 and 7^ =73, 
and we ignore McClure's parameters 7q and 72 .) We 
then have 



H = 



T] 0\ f A B\ fr]"" 
1 



^0 ly \5* A 

with T] = e*'^^i+^2)/3 and 

MOi,0„ 0s) = 73 e-''^ T{0„0,) + 73 e^'^ T*{e„0,) 
B{e„0M = -loT{0M+l^e-^'^T*{0,,e^) 

where 

T(i9i,6>2) = 1 + e^^i + e^^2 
The energy eigenvalues are clearly 

E^{e) = A{d)±\B{e)\ . 

Under a 60° rotation, we have 



(10) 



^2 5 



(11) 

(12) 
(13) 



One then finds A{6') = A(6>) and B(6>') = e''^^ B{6). 
Hence, E^{e') = E^{d). 

Degeneracies identified with a one-parameter family of 
Dirac points occur when B{6) = 0. Solving, we obtain 
the relation 



along the degeneracy curve, where 

7o 7i + 72 73 



Bo 



7i 73 + 7o 72 



= -0.124 



-1.30 X 10" 



The energy along this Dirac curve is 

E{6^) = 8o + W cos (6'i + 6*2 - 36*3) 

with 

£^ = 2r^ 73 = 62 meV 
W = 2r^ 73 = 6.5 meV. 



(14) 

(15) 
(16) 

(17) 

(18) 
(19) 




(^3 = 27r) 



electrons 



{O3 = o) 



FIG. 2: McClure's "sausage link" Fermi surface for rhombo- 
hedral graphite, greatly exaggerated. See also fig. 2 of ref. 

0. 



Since and 7^2 are small, the Dirac curve, when pro- 
jected into the basal Brillouin zone, lies close to the zone 
corners. Note that E{Ojj) goes through three complete 
periods as 6^ advances from to 27r, resulting in Mc- 
Clure's 'sausage link' Fermi surface [10], depicted in fig. 
[21 To find the equation of the Dirac curve, we expand 
about = {6-^,62) = (^7^) at the K point, writing 
6/ = a + C, and find 



T(6)i + (5i9i, 6)2 + (^6>2) = e-^^/^ 66-^ - e 



^^/^(5i92 + 0((56>2) . 

(20) 

Solving for the Dirac line Cl^s) as a formal series in the 
small parameters and r'2, we obtain 



56^ 



2 

2 

V3 



-A sin (193 
sin (6>3 + 



-f)+r2sin (2^3 
f)-r2sin (2^3- 



0{F^) . 



Note that the bandwidth of the Dirac point energies 
is tiny: 2W ^ 13meV. This means that the Landau 
levels are quite narrow - moreso than in Bernal stacked 
graphite. The Fermi surface resembles the sketch in fig. 
[21 which is adapted from fig. 2 of ref. [13] 



A. Weak Fields : Kohn-Luttinger Substitution 

We assume the magnetic field B is directed along z. 
To obtain the Landau levels, we expand about the Dirac 
points. (This is essentially equivalent to expanding about 
the Fermi energy, since the bandwidth of the Dirac points 
is so tiny.) We write 



k 



(21) 
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where tt = p -\- ^A. With 60 j = {k — K) • a-^ we have 
1 



= - TT^ a 



600 



1 

2h 



V3 



(22) 
(23) 



Recah [TT^^TTy] = —ih'^/i'^ where £^ = y^hc/eB is the 

magnetic length. From eqn. [20l to lowest order in 66^ 2? 
we have 



[T, T"^] =2zsin(7r/3) [50^, 
— 27rV^p/q . 



(24) 



where the flux per unit cell area i? = ^ is assumed 
to be a rational multiple p/q of the Dirac flux quantum 
= hc/e. This means we may write 



eb , 



where 



VB/Bo 



(25) 



(26) 



and is a Landau level raising operator: [b^b^] = 1. 
Recah that the field scale Bq = {hc/e)/?>7ia^ = 7275 T. It 
is convenient to define 0^ = ^3 — 1(^1+^2)5 ^ind to absorb 
a phase into the definition of 6, taking T = — ee~*^^ b^ . 
Note that when the magnetic field lies along the c-axis, 
it is exp(zl93) and not ex.p{i6^) which commutes with the 
magnetic translations t{ai 2)- The Hamiltonian is then 



-73(6 + 6+) j^e-'^n"^ -736 
-73(6 + 6t) 

Consider the matrix operators 

e-*^"3 



2o = 7o 



6 + 6t 



2i = 73 I 5t 

The eigenvectors of Qq are 

1*0) 



and 



^ n / 



y/2\±\n- 



1 



(27) 
(28) 

(29) 
(30) 

(31) 



i?o = ±v^e7o, (32) 



where n = 1, 2, 3, . . . . It is easy to see that 

(*^|Qi|*i^) = o, 



(33) 



as well as (^q | Qi|^o) — hence there is no first order 
shift of the eigenvalues. Therefore, up to first order in e, 
the Landau level energies are given by 



EM = £0 + W^cos(3^3) ± , 



(34) 



where n = 0,l,2,... . The gap between Landau levels n 
and n + 1 collapses when 



(35) 



(36) 



e 7o + W = e 7o Vn + 1 - W , 
which gives a critical field of 

^c,n = (Vn+1 + 5c,0 , 

with B^ Q = (2W/7o)2 • ^0 = 0.123 T. 

B. Comparison with Bernal Stacking 



The ABAB stacking pattern of Bernal hexagonal 
graphite is shown in fig. [3l To obtain the critical fields 
in BHG, it suffices to consider a simple nearest-neighbor 
model [5]. Expanding about the K-H spine in the Bril- 
louin zone, we obtain in the presence of a uniform c-axis 
magnetic field. 



H 



I e 7o 6 2 7^ cos ^3 
et, 6t 



V 



2 7]^ cos 6^3 













e7o6 






-e-i^b^ 




(37) 



where e = (2i^^/Zpl qfl^ = y^B/Bo as in the rhom- 
bohedral case. The spectrum has explicit particle-hole 
symmetry. For n = there are eigenvalues at ± (e^ 7q + 




FIG. 3: Top- view of Bernal hexagonal graphite. 
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Magnetic field (T) 

FIG. 4: Landau level structure in rhomb ohedral graphite 
within the tight binding model of section |lll with Zeeman 
term ignored. Principal band gaps are labeled by the Chern 
number C (per spin degree of freedom). When lies within 

a gap, the Hall conductivity is 2C x ^/ (Sd), where d = 3.37 
A is the interplane spacing. 
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Magnetic field (T) 

FIG. 5: Landau level structure in Bernal graphite within the 
nearest- neighbor hopping model, with Zeeman term ignored. 
Principal band gaps are labeled by the Chern number C (per 
spin degree of freedom). When lies within a gap, the 
Hall conductivity is 2C x ^/(2d). When further neighbor 
hoppings are included, particle-hole symmetry is broken, and 
a finite field is required to open the principal gap. 



47^ cos^^g) and a doubly degenerate level at = 0. 
For n 7^ 0, 

^n=(^+i)e' 70 +27^08^^3 (38) 

± \/i 7o' + 4 (n + I) e2 7^ 7^ 008^^3 + 4 7^ cos^O^ . 

In fig. O we plot the lowest several energy bands ver- 
sus magnetic field for BHG. Due to the inadequacies of 
the nearest neighbor model, the principal gap surround- 
ing central E = Landau levels opens immediately for 
nonzero B. Including more realistic band structure ef- 
fects, consistent with the semimetallic nature of BHG, 
this gap opens at a critical field of ^ 15 T for positive 
energies and B^ ^ 7T for negative energies (sj] . The Hall 
conductance is quantized at a^y = 2Ce^/hcQ when the 
Fermi level lies in a bulk gap, where Cq = 2d in BHG 
and Cq = ?>d in RG, where d = 3.37 A is the spacing be- 
tween planes, and C is a topological integer associated 
with the gap. In both cases, the values of C are such 
that CF^y corresponds to the graphene quantization per 
layer, changing by Ae^ /hd as one crosses a Landau level. 
We indicate the width of the bands by shading the re- 
gion between cos^6>q = and cos^6>3 = 1. In both cases, 
the Zeeman coupling is omitted; with g ^ 2 the Zeeman 
splitting is small compared with the cyclotron energy. 



III. CHIRAL SURFACE STATES 

As shown by Hatsugai the Chern number C can 
also be computed by following the spectral flow in a sys- 
tem with edges, wrapped around a cylinder, as a func- 
tion of the gauge flux through the cylinder. To elicit this 
spectral flow, we derive a Hofstadter Hamiltonian [l2[ 
for RG. We start with the Hamiltonian elements in sec- 
tion [III but now treating them as magnetic translations, 
which satisfy the algebra 

t(a) t{h) = e^-xb-/24 t(a + b) , (39) 

where B = 5 n. For our problem we deflne the elemen- 
tary translations 

= t(ffi) , ^2 = t(ai + ffi) (40) 

as well as r = t{dz) = t{a^ + fR). 

Since with n = z we have that r commutes with and 

and we can specify its eigenvalue as e*^3. As for t-^ 2? 
we have 

t,t^=e'^/H^t, , (41) 

where (j) = = 2Tip/q is the flux per graphene 

hexagon in units of hc/e. We may then write 

H=\ , (42) 
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FIG. 6: Spectral flow in rhomb ohedral graphite showing edge 
state evolution. Top panel: armchair edge, perpendicular to 
fR; bottom panel: zigzag edge, perpendicular to ai. The bulk 
gaps are labeled by Chern numbers C which correspond to 
the number of edge states crossing the gap as the angle Oi is 
varied. The flux per unit cell here is rather large, with p = 1 
and q = 200, corresponding to a fleld of 5 = 396 T. The 
topological features of the edge state spectral flow are robust 
with respect to fleld. 



edge, perpendicular to the vector a^^; ^his is shown in the 
bottom panel of fig. [6l For periodic systems, exact di- 
agonalizations performed using the Lanczos method for 
q up to 1500 with the package ARPACK were found to 
agree with the weak field results of section III A[ 



IV. STACKING FAULTS IN BERNAL 
HEXAGONAL GRAPHITE 

We now turn to an analysis of simple stacking faults 
in BHG, first with B = and then for finite B. Consider 
first a triangular lattice, which is tripartite, and its three 
triangular sublattices A, and C. By eliminating one 
of these three sublattices, the remaining structure will be 
a honeycomb lattice. Now imagine a stack of triangular 
lattices. At each layer, we choose a sublattice A, 5, or 
C to remove; this defines a stacking pattern. Since it 
is energetically unfavorable to stack a honeycomb layer 
directly atop another, at each layer we have two choices 
consistent with the layer below. If the empty sublattices 
are in ABC et eye. order from layer / to layer / + 1, 
we write 0-^1 n-\-i ~ -"-^ instead the order is CBA et 
eye.^ we write cr^ = —1. For RG, the a indices are 

'ferromagnetic', i.e. + + ++ or . For BHG, the 

indices are ordered 'antiferromagnetically', i.e. H \ — . 

The three three triangular sublattices A, B, and C are 
defined by 



with 



<,,n2 = ^1^1 +n2a2 +ffii 



(47) 
(48) 
(49) 



73 (e^'^ ti+e 



(43) We define three additional sublattices by 



-i(/)/6 



and 





^2 




"2 




77,2 


= ni ai - 


^ ^2 ^2 - 


hffii 


4) 


4 = ^BB 




"2 


= < 


77,2 


= ni ai - 


^ ^2 ^2 - 










"2 




^2 


= ni ai - 


Vn^ a2 





H 



-2103 



AB 



ill e*^3+72e 

+ (736-''^ 
-(70 

We define the basis {| ^)} as follows: 



7o 4 + 73 e 

i(f)/6 



7oe 
73 e'^' e-^^/^ 



t 



t^\n) 



gm0/3 ^ 



) 



I n) = e''^^ I n + 1) , 



(44) 

(45) 
(46) 



(52) 



The sites {'?^a(^i? ^2) ? '^a(^15 ^2)} ^^c- fo™ a honey- 
comb lattice, which we call the A or a structure. Bernal 
graphite is stacked in an ABABAB configuration. 

Within each honeycomb layer, we write the wavefunc- 
tion as a two-component spinor, 



where 0-^ = O-^/Sq and | n + 3^) = | n). Taking the ma- 
trix elements of H within this basis, one obtains a rank 6q 
matrix H to diagonalize, with periodic boundary condi- 
tions. If we introduce an edge by eliminating the coupling 
between states | l) and | 3^), and plot the spectral flow 

as a function of O-^^ we obtain the top panel in fig. [6l We 
can also obtain the chiral surface state flow for a zigzag 



(53) 



where k is the crystal momentum in the basal {kz — 0) 
Brillouin zone. 

The hopping between planes is described by the follow- 
ing local Schrodinger equation, which couples a central 
plane / to planes below (/ — 1) and above (/ + 1): 



G and 



g' ^ i.e. the shift in the u 
sublattice sites from plane / — 1 to plane / is through a 



Here, cr^_i , 
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vector crfR^. The matrix M is given by 



From /C ^ = 0, we may write 



and 



and 



M-( ^ 



1 






1 0/ • 



(55) 



(56) 



(57) 



A. Bernal Hexagonal Graphite 

We first consider the BHG stacking order ABABAB^ 
where cr nj^i = (—1)^- Using translational invariance, we 
may write, for the even and odd sites 



2j 



2j+l 



where 



M(/) + cos {k^d)E-X = 

MX + cos {k^d)E^(l) = . 
Inverting the second of these equations gives 

X = -27i cos {k^d)M-^U^(l) . 
Substituting this into the first equation yields 

Accordingly, we define 

JC = M -4jlcos^{k^d)E-M-^E^ 
/ E 7o^ 



(58) 
(59) 

(60) 
(61) 

(62) 
(63) 

(64) 
(65) 



Setting det/C = yields the eigenvalue equation for 
Bernal graphite, 

{E^ - 7o' \S\'f - AE^ 7? cos2 {kj) = , (66) 
with solutions 



e: 



= -M7i cos {Kd) - MV7?cos2(fc,d)+7o'|5|2 , 

(67) 

where /i = ±1 and ji' = ±1. The four choices for (/i,/i') 
correspond to the four energy bands. 



From eqn. [62| 



then, we have 
2£^7;l cos{k^d) 



E 



-E 



(68) 



E 



(69) 



B. Step Defect 



Consider now the stacking defect ABABCBCB, which 
in terms of the variables may be depicted as 

■■•l + l-l + l-l-l + l-l + l--- (70) 

The central plane we label / = 0. For plane indices / < 0, 
the odd layers correpond to (j) planes and the even layers 
to X planes. For / > 0, the even layers correspond to (j) 
planes and the odd layers to X planes. With / < 0, we 
consider an incident plane wave running to the right (up) 
and a reflected plane wave running to the left (down). 
Then we have 

V^2i = e^'^'^^"^ + f^' e-2^^'^^^) X (71) 
V^2,+i = e^^'^'+'^'^' + f3' e-^(2^-+i)^^^) (/) , (72) 

for all j < 0. Here a is the complex amplitude of the 
incident wave and (3' is the complex amplitude of the 
reflected wave. 

Correspondingly, we have 



(^g2ii/c,d^^/g-2ii/c,d>j^ ^ (74) 



for all j > 0. Here a' is the incident amplitude (from the 
right/top) and P is the reflected amplitude. 

To match the solutions for positive and negative /, we 
first invoke eqn. [SHwith / = — 1: 



The most general solution for i/jq is then 
^o = (a + /3')X+(^j) , 



(75) 



(76) 



where b is an arbitrary complex number. Note that U~ 
annihilates any vector with upper component 0. 
Next, set / = +1 and obtain 



MV;i+7i^+V^o + 7i^+V^2 
We may now write 







% = (/3 + aO(/> + 



(77) 



(78) 
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where a is an arbitrary complex parameter. Note that 
annihilates any vector with lower component 0. 
The parameters a and b are then fixed by equating 
these two expressions for ?/^q, yielding 



(^5) =(" + /?') X- + 



(79) 



The wavefunction at / = can now be found. One simple 
way is to take the upper component from eqn. 11611 and 
the lower component from eqn. 11621 



(80) 



Next, we write the Schrodinger equation for the I = 
plane: 



_fE ^oS\ ({a + (3')x; 



(81) 



7o5* E 

+ 7i (ae-''=^^ + /3'e 

+ 7i {13 e'^'" 







This yields two equations which may be solved to relate 
the outgoing amplitudes /3 and /3' to the incoming am- 
plitudes a and a\ i.e. to derive the 5-matrix. Using our 
previously derived results for cj) and X, we find that the 
above equation reduces to 

"-[joS* E [{P + a') 



+ 7i 



(ae" 
,M {13 i 



-ik^d 
ik^d 



■[5' e 



a e 



I ^ik^d\ 



-ik^d\ 



(82) 



This yields 




/ 7o5 fiE + j,e^o^/^ 

where 0^ = 2k^d. The <S-matrix is defined by 

*S— matrix 



t r' 
r t' 



(83) 



(84) 



Solving for S, we obtain 

_ ^2i7o5*sin(e,/2) 

S 



7i 



7i 

2i7o5sin(e,/2) 



7i cos(0,) + 2iii sm{e,/2) 7^ cos2(0,/2) + 72 |5|2 

(85) 
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FIG. 7: Reflection and transmission coefficients for k 
cx\ bi + a2 b2, for four sets of (ai, 0:2), in the vicinity of the 

^ 3 ' 3 / 



Dirac point (|, |). Only positive energies are shown. 



Thus, for all /i and /i', we have 



and 



R 



T 



7i 



7? + 47o'|5Psin2(^,/2) 

47g|^psin^(^,/2) 
7? +47o'|5Psin^(^,/2) 



(86) 



(87) 



As k approaches either zone corner K or K\ the trans- 
mission goes to zero. This is because the chains which 
extend through BHG are cut and shifted at the stacking 
fault. Curiously, the transmission coefficient T goes to 
unity when 7^ ^ 0. Note also that along K-H and K'- 
H' we have 5 = and hence = 1, T = 0. At the band 
edges, we have 



R{e, = {)) = ! , R{e, = ^) 



7? 



(88) 



with T = 1 — for the transmission coefficients. 

C. Existence of Bound States 

To search for bound states, we take, for j > 0, 

V'„=_2, = X V„=2, = 13 <t> (89) 

V„=-2i+i=e™</' ^„=2i+i=/3e-™X, (90) 

and solve for and E. At the plane / = we have 



9 



0.25 




H 0.1 



0.05 [- 




f////////////A 



_L 



0.02 



0.04 0.06 



0.08 



0.1 



FIG. 8: Energy bands and bound state dispersion for 71 = 
0.1 70 for small values of \S\. The bulk bands E^'^^ are de- 
picted by the red and blue hatched regions, respectively. The 
bound state is the thick black dot-dash curve. 



The Schrodinger equation for / 7^ then yields 

M X + 271 cosh(A>:) ^7+0 = (92) 
M ^ + 27i cosh(/^) i:"X = . (93) 

This yields 



E = -/i7i cosh(A>:) - ^' ^-fl cosh'(/^) + 72 |^|2 , (94) 
where once again /i = ±1 and 11' = ±1. Again we have 



E 



At ^ = we again have 

MVo + 7i^+V'-i+7i^"^+i = , (96) 
which here yields 



(97) 



This yields two equations for which may be written as 



/3 



(98) 



This fixes the energy at 

E = -^,^^e-^±^^\S\. (99) 

Thus, we have a bound state at positive energy (and 
a corresponding one at negative energy) for each real, 



positive value of z^, which solves one of the four equations 
(for /i, fi' = ±1) 

-IJL^i e""" + /i'7Q = -/i7i cosh(A^) 



-MV7?cosh2(«)+7o'I^P • 



(100) 

We assume 7q > 0. In the SWMC analysis [3], one has 
7i ^ — 390meV and 70 ~ 3.2 eV. The vertical hopping 
is negative due to the sign of the overlap of orbitals 
on consecutive layers. In order to have a bound state 
solution, we must have fifi' = sgn(7;L)5 resulting in 



V't? cosh2(«) + li \S\^ - 7o \S\ = |7i I sinh(K) , (101) 
the solution of which is 



sinh(A>:) = 



7i 



2 7ol^l 



: u . 



(102) 



Thus, there are two bound states for all k in the Bril- 
louin zone, one at positive energy, corresponding to the 
choices ji = ji' = sgn(7;|^), and one at negative energy, 
corresponding to the choices ji = ji' = — sgn(7]^). Solving 
for we have 



e'^'^ = ±u^^l^v? . (103) 
The bound state energy may now be written as 

1 



|7i| + 



2u 



+ (9(ii-^) 



(104) 
(105) 



where the expansion in the second line is for large i.e. 
^q\S\ <C |7i|. Note that the bound state disperses as 
|kp. Recall for Bernal graphite that the dispersion is 
linear in |k| in the vicinity of H and quadratic elsewhere 
along the K-H spine. The length scale associated with 
the bound state is . For ^ oo, ~ l/ln(2^i). 

Since the spectrum, including bound states, is particle- 
hole symmetric, we may without loss of generality limit 
our attention to £^ > states. The continuum bands, for 
fixed k, range over energies 



7i 



7ol^l <^r<|7ih 



(106) 



It 



(107) 



The bound state we have analyzed lives just below the 
bottom of the E^^ band. The binding energy is = 



n(3) 



l7il 



£^3, and is given by 
1 



+ 4^/2 - l) + VTT^ -l-u . (108) 
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FIG. 9: Binding energy of the bound state versus \n(u), where 
u = |7i/27o6'|. 



In figs. [8] we plot the bound state spectrum for the case 
ITiI/To — small values of 1*5 1, i.e. close to the zone 

corners, where u is large. At the zone center, IS*] = 3 is 
maximized and u achieves its minimum value; for refer- 
ence, i^gwMc — 0.02057. The binding energy vanishes 
in both the u ^ and u ^ oo limits, as shown in 
fig. [9l The maximum of A occurs for = 1, where 
^/ITiI — 0.03225, corresponding to a binding energy of 
approximately 13meV. In the appendix, we compute the 
bound state spectrum for the full SWMC model. 



V. FINITE B CASE 

To obtain the Landau levels, we expand about the 
Dirac points, following the method described in section 
III A[ We have ^K+fi/fi ~ ^^^^ ^ given in eqn. [26l 
At 5 = lOT one has e = 0.0371. With 7^ = 0.39 eV 
and 7o = 3.16 eV, we have r = 0.123 and e/r = ^J B jB-^ 
where B^ = 110.8 T. For physical fields, then, we have 
6 < r. Note that one can also write 



^70 



(109) 



where = -/Z^QajlK is the Fermi velocity (a = 2.46 A 
is the lattice spacing in the hexagonal planes) and 1^ = 
hc/eB is the magnetic length. 



A. Bernal Stacking and Landau Levels 

We define the operator- valued matrix 



M 



E 6 
£70 6t E 



For perfect Bernal stacking, we have 

Mi^^j + 7i ^+(^2i-i + ^2i+i) = 
M4'2j+i + 7i {^2j + ^2^.+2) = 



(110) 



(111) 

(112) 



We now write the wavefunction in terms of right and left 
moving components: 



iP^j = (/e'« +0'e- 



a I n) 
P\n+1) 



(113) 



(114) 



where we assume n > 0. We therefore have 

M„ Q + 27i cos(g/2) S- (^^^ = (115) 
M„+i (^^ + 27i cos{qj2) E+ Q = (116) 



where 



/n e7Q E 



(117) 



This leads to 



P„(i?) = det M„+i - 47? cos'{q/2) S+ M"! 



A-iiE^ cos\q/2) 



E"^ - ne^ 7o 



• (118) 



Setting Pn{E) = yields the spectrum E = E^{q) of 
Bernal hexagonal graphite: 



^n,±(^) / , 1\ 2 I o 2 2/ /o^ 

— '-^ — = (n + ^) e + 2r cos [q/2) 
7o 



(119) 



± ^/4r4 cos4(9/2) + (4n + 2) r"^ cos^{q/2) + , 



where r = 7i/7o- Expanding for small e/r, we have 



\/n{n + 1) — , fiVn 



Vn+T € , 2r + (n + i) — + . . . 



(120) 
(121) 
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Magnetic field (T) 

FIG. 10: Landau levels in graphite. Subbands £^n=i,+ (red), 
En=i,- (blue), and En^o,-\- (green) are shown. The zero 
modes are shown in black. 



B. Zero Modes 

The case n = must be considered separately. Con- 
sider the wavefunction 



2i 







r , ■>p2j + l 



(122) 



This is an £^ = eigenstate for any J. It describes a 
state localized on a single plane. 

We can find more solutions by writing 



= 



The Schrodinger equation then requires 



ve7o E 
on even planes and 



= 



e(0)+27,cos(,/2) (0 



(123) 
(124) 

(125) 

(126) 



on odd planes. Thus, we have three equations for the 
remaining three eigenvalues: 

= £; a + e 7o /3 + cos(g/2) y (127) 
{) = E[5^e-fQa (128) 
= ^ y + 27i cos((7/2) a . (129) 
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FIG. 11: Landau level indices for scattering at a stacking fault 
in Bernal graphite. 



We immediately see that £^ = is an eigenvalue, with 
eigenvector 



/ 

/3 = -27iCos(^/2) 



(130) 



If we Fourier transform this solution, multiplying by 
g-*q(j+2) Q^T^^ summing over we find a purely local- 
ized state, with 



V^2J = -7i 
V^2J+l =^7o 



2J+2 



-7i 






0). 


ii>; ' 



(131) 
(132) 
(133) 



with all other t/j^ = 0. This zero mode is localized 
on three layers. The remaining two solutions are easily 
found to be 



E 

-e7o 
-27i cos(g/2)^ 



(134) 



with E = Eq^^ = ±\/e^7o +47i2cos2(g/2). These solu- 
tions are wave-like and disperse with q. 



C. Stacking Fault 

For the system with a single step stacking fault, the 
situation is as depicted in fig. [TTl We then swap the 
notation for even and odd planes on the right half of 
the system (layer indices / > 0) with respect to eqn. I114[ 
and introduce wavevectors q^ and q^ for the left and right 
half-systems. We must then match the energies on left 
and right sides of the fault: 



(135) 
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To identify the bound states, we write the wavefunc- 
tion for / > as 



a - I n + l) 



"^2^-1- U.|n + 2) 



and for / < as 

/ - 

At / = we write 



2i 



I n) 



(136) 



2i 



I "> 



Xo\n) 
yo\n+l)J 



(137) 



(138) 



The Schrodinger equation, evaluated for both even 
and odd planes with / > and / < now gives 
eight relations among the eight sets of coefficients 
{ a J , Pj , Xj , yj , a J , Pj , Xj , } , expressible as 



E 



VnTT E 



+7i 







and 



VnT2 670 ^ 



-7, 



and 



= 
(139) 

= 
(140) 



E 



Vn + 1 67o\ 



VnTT 6 7o 
and 



-7i 



(141) 



E 67o 
/n 6 7o £^ 



■7i 



. 



(142) 

We can use these equations to eliminate the four sets 
of coefficients { /3j , Xj , ^Vj}' 



^3 



^3 



= -Vn + 2 6 7o a. 
yj = -VnTT 6 7o Xj 
-VnTT e-fQE-^y. 

We then obtain 

= ^n+l(^)?/j +«j 

= ^n+2(^) (^j^Vj-l^Vj 

= i?,(^)^^. , 



where 



(143) 
(144) 
(145) 
(146) 

(147) 
(148) 
(149) 
(150) 

(151) 



with = n6^ 7o. We then have 



a 



^0 



and 



/5i 



where 



KniE) 



^ fRMRn+iiE)-l R„{E) 

-Rn+liE) -1 



(152) 



(153) 



(154) 



Note that deti^^(£^) = 1, and that the characteristic 
polynomial det (A — K^) is real for real A. It is easy 
to see that the eigenvalues of K^{E) form a complex 
conjugate pair e^*^ if the energy E satisfies the condition 
the condition |Tri^^(£;)| < 2, or 



0<R^{E)R^^,{E)<4. 



(155) 



This is the condition that E lies within one of four energy 
bands. The roots of R^^iE) Rn+i{E) = lie at E'^ = El 
and E'^ = E^^-^, while the roots of R^iE) R^^^iiE) = 4 
lie at E'^ = _ and E'^ = E^^^^^, where 



<± = (n + ^)e'7o+27f 



(156) 



± ^4 7f + (4n + 2) £2 72 ^2 + 1 ,4 ^4 _ 
The bands are then given by 

El_<E^<El , El+,<E^<El^. (157) 
In the limit a = 6^ 70 /71 <C 1 , we can expand and write 



n(n + 1) 



4 4 



47? 

~472 + (2n+l)e272 . 
At plane n = the Schrodinger equation yields 

XoJ V7i E J 



(158) 
(159) 



7i 

a/tT+T e7o 



(160) 



D. Scattering Matrix 



If both |Tri^^(£;)| < 2 and \Tt K^^^{E)\ < 2, then 
we can write 



^1 



and 



1 1 ^ +o*^""^^\ 



(161) 



(162) 
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FIG. 12: Bulk energy bands (shaded and hatched regions) and 
bound states (magenta curves) versus magnetic field for 70 = 
3.16 eV and 71 = — 0.39eV (tight binding; nearest neighbor 
hopping only). The lowest energy bound state merges into 
the band continuum at B ^ 15 T. The other bound states 
remain sharp over the energy range shown and do not mix 
with the lowest bulk band. 



where 



(n) 



(163) 



Then we have 



(164) 



(165) 

The tS-matrix, which relates incoming to outgoing flux 
amplitudes, is then obtained from eqns. 1160^ I16H and 
11621 upon replacing / v\!^X^ O 0'v\!^ ^ I' 
v\l^X^ and O Ov\l^ ^ where = dEnici^jdq^ and 

E. Bound States 

If a state is evanescent on both sides of the stack- 
ing fault, we must have that both |TrK^(£^)| > 2 and 

\TrK^^-^{E)\ > 2. The eigenvalues of K^{E) are given 
by 



A, 



n,± 



Ir ± -\/t^ - 4 

2 'n ^ 2 V 'n ^ 



(166) 



where 

rJiS) = Tr ifJiS) = R^(E)R^^^{E) - 2 . (167) 



In order that the solution in eqn. 11531 not blow up for 
n ±00, we must require that \ and \_} \ have 

no weight in the > 1 eigenspaces for K^{E) and 
K^_^-^{E), respectively. This means 



2+1 Vo 



(168) 
(169) 



where < 1. When we combine these equations with 
those in eqn. I16U1 we obtain 



M 



Ji, 



(170) 



where 



M 



7i 


V 



E 







n+l 







\^A< 
E 



I 



(171) 



Vn + 1 7i 

A solution requires D{E) = det A^(£^) = 0. We have 

Z)(i?) = [£;i?„+2-7i(l+^<+i)] [i?i?„+i-7i(l + ^;?)] 
-(n + l)e2 72i?„^^i?„^2 . (172) 

Let us look for a bound state with energy E which 
is parametrically (in a) smaller than both and e7o. 
Then Rn{E) ~ -E^/Ej-^^, from which we obtain yl< ~ 

(i?„i?„+i)"'- Then find 



D{E)^-/f-in + lfin + 2) 



£676 



7fE2 • 

Setting D{E) = yields the bound state energy, 



E' 



{n+lf{n + 2) 



(173) 



(174) 



Thus, the bound state energy is proportional to 

In fig. [T2I we plot the lowest ten bound state energies 

versus magnetic field. 



VI. SURFACE SPECTROSCOPY OF BURIED 
STACKING FAULTS 

Our previous results for the transmission through a 
stacking defect suggest that these defects are very effec- 
tive in decoupling graphene stacks. We analyze now the 
density of states at a graphite surface in the presence of 
a stacking defect a few layers below the surface. The 
stacking sequence is {KB) ^^i^^^C^- • • . The number of 
layers between the surface and the defect is N. 
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FIG. 13: (Color online). Left: Density of states for the two 
inequivalent sites of a graphite surface with a stacking defect 
20 layers below the surface Triangles (red) give the density of 
states at the site with a nearest neighbor in the layer below. 
Squares (green) give the density of states at the site without 
nearest neighbors in the layer below. Right: as in the left 
panel, with a defect 100 layers below the surface. 



The system can be separated into a perfect semi- 
infinite graphite sample coupled to the defect layer, and 
N layers between the defect and the surface. We will 
only include the parameters and . The semi-infinite 
portion can be integrated out. The site of the defect layer 
connected to it acquires a self energy: 



(175) 



We now integrate out this site, leading to the self energy: 



uj - Eo{uj) 



(176) 



The procedure can be iterated leading to new self energies 
for sites 2, 3 . . .N, resulting in the hierarchy 



\loS\' 



7? 



(177) 
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FIG. 14: (Color online). Surface density of states for a semi- 
infinite stack with a defect ten layers below the surface. The 
Landau level index is n = 2, and the fields studied are B = 1 
T (red) and B = 10 T (blue). Left: Sublattice with a nearest 
neighbor in the contiguous layer. Right: Sublattice without 
a neighbor in the contiguous layer. 



surface layer {N) are: 

C(^) = - 



and 



(178) 



(179) 



We show in fig. [13] the surface density of states when 
such a defect lies twenty and hundred layers below the 
surface, obtained by integrating the imaginary part of the 
Green's functions in eq. 1 1791 over the in-plane component 
k|| of the wave vector. 

The density of states show a number of resonances, 
which are smoothed out when the number of layers be- 
tween the defect and the surface is large. For N ^ 1^ we 
recover the analytical results in ref. [ij]. These results 
are consistent with the analysis in the previous sections, 
which show that the transmission through the defect is 
strongly suppressed. The layers between the defect and 
the surface become effectively decoupled from the bulk 
of the system. 

The previous analysis can be extended to the study 
of Landau levels in a magnetic field. As discussed ear- 
lier, the hoppings within the layers depend now on the 
Landau level index n, instead of on ky. The n depen- 
dence of the hoppings in he two layers within the unit 
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FIG. 15: Surface density of states at the sublattice without a 
nearest neighbor in the next layer. The system has a stacking 
fault of the type described in the text ten layers from the 
surface. Top: n = 2. Bottom: n — 10. 



cell is different. Because of this, the self energy obtained 
by integrating out the perfect semiinfinite region leads 
to a more complicated expression than those in eq. 11751 
Within the region between the defect and the surface the 
successive self energies have a twofold periodicity: 



7? 



-2 
B 



{n-\)vl 



7? 



(180) 
(181) 



The resulting densities of states for Landau level index 
n = 2 and two magnetic fields, 5 = 1 T and 5 = 10 T, 
are shown in fig. [TH 

We show finally in fig. [15] the dependence of the peaks 
in the surface density of states on the magnetic field. As 
before, there is a stacking defects ten layers below the 
surface. In agreement with experiments [H, HB, [13 ^ there 
are peaks which scale as \fB and peaks which scale as B. 
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FIG. 16: Sketch of the expected behavior of the Hall conduc- 
tivity in the quantum limit, for a lightly doped system. The 
leftmost plateau is a bulk effect, related to the Landau levels 
of Bernal graphite @|. The size of the other jumps depends on 
the concentration x of stacking defects. The continuum bands 
of Landau levels lead to a monotonically varying conductivity. 



VII. DISCUSSION 

We have analyzed the appearance of two dimensional 
features in bulk graphite. We show that deviations 
from the Bernal stacking order are very effective in in- 
ducing two dimensional behavior. An ordered array of 
graphene layers with the rhombohedral stacking order 
leads to isolated Landau levels, and to quantized quan- 
tum Hall plateaus at moderate magnetic fields in doped 
systems. We found that the gap between Landau level 
subbands of indices n and n + 1 opens at a field 5c, n with 
5c,n=o = ^0 = 0.123 T and 5c, n ^ ^nB^ for large n. By 
contrast, in Bernal graphite, the first gap is predicted to 
open at fields on the order of 10 T [5], and the second gap 
opens only at enormous field, on the order of 1000 T. 

We have also considered the simplest stacking defect 
in Bernal graphite, which has locally a rhombohedral ar- 
rangement. These defects are expected to be common 
in many graphite samples, and concentrations up to 10% 
have been reported [H . These defects are very effective in 
decoupling the electronic states on either side. They also 
give rise to a two dimensional band of electronic states, 
localized in the vicinity of the defect. Within a nearest 
neighbor tight binding model for the tt band of graphite, 
with in-plane hopping 7q = 3.16 eV and interplane hop- 
ping 7^ = 0.39 eV, we found a maximum binding energy 
of approximately 13 meV, for states rather close to the 
corners in the basal Brillouin zone. When the full SWMC 
model is taken into account Q, we obtain a maximum 
binding energy of almost 40 meV for electron states and 
20 meV for hole states; the binding energy is significant 
only along the zone faces. 

What are the implications of our work for magneto- 
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transport in graphite with stacking faults? To describe 
the physics, it is helpful to keep in mind the bound state 
Landau level structure of fig. [121 First, suppose the 
graphite is undoped. In this case, the Fermi level remains 
pinned within the central n = Landau levels. With only 
nearest neighbor hoppings, there are two flat bands {i.e. 
which do not disperse as a function of k^) at £^ = asso- 
ciated with each zone corner in the basal Brillouin zone. 
Taking into account the weak second-neighbor plane hop- 
pings 72 and 75 , these bands disperse and acquire a width 
of about 40 meV. For the full SWMC model, due to the 
breaking of electron-hole symmetry, the Fermi level can 
drift within these central Landau subbands, even if the 
system is at electroneutrality. As shown by Yoshioka 
and Fukuyama [23], due to interaction effects one then 
expects a charge density wave (CDW) at sufficiently high 
fields. Anomalies in the observed magnetotransport data 
corresponding to this CDW transition have indeed been 
observed [21]. The presence of stacking faults, which pro- 
duce bound states away from the central Landau levels, 
should not affect this picture. 

However, if the graphite is lightly doped, a different 
picture emerges [5]. In this case, the central Landau lev- 
els become filled at a field B* = ^nd(j)Q^ where n is the 
bulk carrier density, d is the interplane separation {i.e. 
the c-axis lattice constant is 2d due to Bernal stacking). 
For fields B < B"" ^ the central Landau bands are filled, 
and the Hall conductivity should be quantized at a value 
2e^ /hd |5|]. As B is decreased further, the Fermi level 
crosses the bound state energy. The bound state Lan- 
dau levels (one for each spin value and inequivalent zone 
corner) then makes a contribution to a^^, of magnitude 
Acrjf^^* = 2xe^ /hd^ as shown in the sketch in fig. [161 
where x is the concentration of stacking faults. Upon 
further reducing 5, the Fermi level enters into the first 
bulk band, and a^y begins to rise continuously. As 
crosses other bound state Landau levels, additional small 
jumps of Aa^^y^^ = 2xe^ /hd should appear. At a finite 
concentration x of stacking faults, the bound states will 
themselves form a band, and the small jumps will no 
longer have infinite slope. 

The scenario discussed here shows how anomalous fea- 
tures could occur in the high field magnetotransport of 
doped graphite, however we cannot find any obvious con- 
nection between our work and the observations of Kempa 
et al. [H. 

Stacking defects below a graphite surface decouple 
the surface region from the bulk, leading to quasi-two- 
dimensional behavior, with localized Landau levels. We 
have shown how such buried defects leave a signature 
which can be measured by surface spectroscopy. 

Finally, our results suggest that the electronic proper- 
ties of few layer graphene samples can be substantially 
modified by changes in the stacking order. 
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IX. APPENDIX : FULL SWMC TREATMENT 
OF STACKING FAULT 



We define the vectors 



{n < 0) 



(n > 0) 



(182) 

For a stacking defect ABABCBCB the SWMC cou- 
plings are depicted in fig. \17\ In fact, additional cou- 
plings must be introduced at the defect. In the bulk, sites 
have either zero or two c-axis neighbors, but at the stack- 
ing fault there are two sites with a single such neighbor. 
One expects the associated on-site energy A'' ^ ^A. In 
addition, there are three interlayer couplings at the de- 
fect which in principle are distinct from 73 and 74, and 
which we denote in the figure by dotted pale blue lines, 
with hopping amplitude 74. For simplicity, we shall take 
A'' = A and 74 = 74 for two of the links, and 74 = 73 
for the other link. For details, see the definition of the F 
matrix below. 

Let each pair of layers be indexed by a nonzero integer 
n. From the figures, we can read off the Schrodinger 
equations 



where, 

IS], 



K 



(n < -1) (183) 
1 ^ v^, ^ ... v^^+i - (n > 1) , (184) 
consistent with the full SWMC Hamiltonian 0, 



and 



/ -E 
\ I3S 



M = 



-loS 
-E + A' 

7i 
74 5* 



I4S 

-E + A' 

-7oS* 



73^*^ 
74 5 

-E 



(185) 



(b2 





74^ 


73 





575 


7i 


74 5 








375 





V 








572 / 



(186) 



We take the SWMC parameters from ref. [18| 

7o 3160 meV 7^ 390 meV 

72 = — 20meV 73 = 315 meV 

74 = 44 meV 75 = 38 meV , 
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FIG. 17: SWMC couplings for a stacking defect in Bernal 
graphite, showing more clearly the four sublattice structure 
on either side of the defect. 



with A = — 8meV. Here, A' is a combination of the 
original SWMC parameters: 



A' = A + 73 - 72 , 

hence A' = 50 meV. 

At the defect, the Schrodinger equation yields 



where 



fb2 73^ 74 5*\ 

74 5* 7i 

175 

V i72 y 



(188) 
(189) 



(190) 



A. Scattering matrix and bound states 

We write tp^ = z"" X (for n < 0) and (j)^ = z*"" X* (for 
n > 0). In the bulk (|n| > 0), we then have (for both 
sides) 

{z-^M + K^zM'^)X = . (191) 

In order for a solution to exist, we require 

P{z) = det {z-^M + + z Mt) = , (192) 

which is an eighth order equation in z. Note that P{z) = 
guarantees that P{z*~^) = 0. It can also be shown, 
due to the form of M, that P{z) = P(^"^). Thus, the 
allowed values of z come in sets (z, z*, z~'^, z*~^). 

Within a bulk energy band, two of the eight z roots are 
unimodular, and may be written as z^ = e*^ and z^ = 
e~*^ with k real. Their associated eigenvectors are Xi 5. 
Of the remaining six roots, three (2^2, z^^ z^) each have 
modulus greater than unity and are thus unnormalizable 



on the right. The remaining three roots (zg, Zj^ z^) each 
have modulus smaller than unity and are unnormalizable 
on the left. We keep only the normalizable solutions and 
write 



n<0 : V^^=Te^^^Xi + 0'e 



(193) 



^2 ^2 -^3 ~^ ^3 ^3 '^3 ~^ ^4 ^4 '^4 



n > 



J', 



^6 ^6 



5 (194) 



Equations (|188ll89p then yield eight equations in the 
ten unknowns {X,0,T ,0') and (A2, A3, A4, Ag, A7, Ag). 
These then determine, for each energy the 5-matrix, 
defined by the relation 



t r' 
r t' 



(195) 



If two bands overlap, then we have eigenvalues 5 = 



(187) 6^*^ ^2,6 — e^*^, with \z^^^\ > 1 and \zy^s\ ^ 1- 



n < 



ih 



O' e-'^" X5 (196) 
2 + O' e-'^'" Xe 



^3 ^3 -^3 + ^4 ^4 



n > 



= J' e"*" Xt + O e*" x; (197) 
+ f e-^f" Xj + O e'P" X^ 

I A ^^'^ V* _L A r/*^ V* 
-j- Ay -p ^3 ^3 A3 . 

The /S-matrix is then 4x4, and we should take care to 
properly define it to act on flux amplitudes^ viz. 



O' 



1/2 



= S 



1/2 
'h,k 



(198) 



where '^1 = dEi{k)/dk and '^2^^ = dE2{p)/dp. If three 
bands overlap, the 6'- matrix is 6 x 6. 



B. Bound states 

When E does not lie within a bulk band, we write 
n < 



^^=A,z^X,^A2Z^X2 (199) 
-|- z^ X3 H~ ^4 ^4 X4 



n > : 



A5<"x;+^6^6*"X; (200) 



% • 

^7 ^7 A7 -p ^3 ^3 A3 . 



Here, l^^^ 2,3,4] > 1 and |2;5 g ^ g] < 1. Without loss of 
generality, we may assume 



z, 



■n+4 5 



(201) 
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for u = 1,2,3,4. Equations (|188ll89p now give eight 
homogeneous equations in the eight unknowns ^i_8- A 
nontrivial solution can only exist when the correspond- 
ing determinant vanishes, which puts a single complex 
condition on the energy E. The solutions are the allowed 
bound states. 

We now apply eqns. 11881 and 11891 



to 



/=5 



usmg 



(z-^M + i^ + zM'^)X = 



This yields 

which, expanded, gives 



1=5 



(202) 
(203) 



(204) 



(205) 
(206) 



=M*(/)o , (207) 



(208) 



Note that det((5) = det(M) • det(Ar) = 1 independent 
of K and of the above-diagonal elements of M and the 
below-diagonal elements of N . From row reduction, it is 
easy to derive 



N 



7275 



1 75 







74^* -7i7275 ^ h2 
V-72"S375 ^ 74^* \-f,) 



\ 






(214) 

The equations (|208p and (|209p can now be written as an 
8x8 system. 



= , (215) 



where a, 6, and u run from 1 to 4, and / runs from 5 to 
8. The implied sums for each submatrix are over h and 
not u OT I, and is the matrix of eigenvectors of Q: 




kj 



1 



(no sum on i) , (216) 
The bound state 



and 



where i, j, and k run from 1 to 
condition is det(i?) = 0. 

We have numerically found the bound states lying in 
the gap between the bonding and antibonding pi bands 
of graphite. Our results are displayed in fig. [181 For 



1=5 



(209) 



These give 8 homogeneous equations in the 8 unknowns 
can only be solved when the corresponding determinant 
vanishes, which is the condition that E lie at a bound 
state energy. 



C. Method of solution 

Eqn. 11911 can be written as two coupled equations, 

z-^MX + X' = (210) 
i^X + zM'^X-X' = . (211) 

These equations may be recast as the rank-8 system, 

'z^NK -N\ fx 



M 



z \X 



, 



(212) 



where N = ^ . Thus the solutions Zj are the complex 
eigenvalues of the matrix 



-NK N\ (X 
-M I \X' ) ' 



(213) 
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FIG. 18: Binding energies for bound states within the gap at 
negative energies (blue short dash - dot curve) and positive 
energies (red long dash - dot curve), for the SWMC model 
with a single stacking fault, as a function of wavevector in 
the basal Brillouin zone. The solid green line shows the en- 
ergy gap between bonding and antibonding tt bands, which 
collapses in the vicinity of the basal zone corner K. 
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the full SWMC calculation, there is no longer particle- 
hole symmetry. We find that the binding energy {i.e. 
the distance of the bound state from the closest band 
extremum) is considerable along the entire KM edge. 
This is in contrast to our analytic results for the nearest 
neighbor model, where the bound state energy was con- 
siderable only for l^*! ~ 7^/270 0.062, which is satisfied 
only on a small ring about the K and K' points. On the 
other hand, the lack of bound states along the FK edge 
implies a finite broadening of the Landau levels derived 
from this band. 

It is important to realize that the SWMC model itself 
is only valid close to the K-H spine in the Brillouin zone. 
The model must be extended, as in ref. [l^, to include 
other tight binding parameters, in order to fit the tt band 
throughout the entire zone, which is necessary in order 
to model various optical transitions. In this case, the 
in-plane hopping is modified: 

7o5^7^''5i+7^''52+7^''53 , (217) 



where 7^^^^ and are, respectively, the amplitude and 
lattice sum of e^^'^ corresponding to the n*^ nearest 
neighbor in-plane inter-sublattice hopping [l^, subject 
to the constraint 

^swMc^^a)_2^(^)+^(3) . (218) 

It is a rather simple matter to include such effects in 
our calculation, and we find in general, for a broad set of 
possible parameterizations satisfying the constraint, that 
our results have the same qualitative features. 

Our approximations regarding the parameters A'' and 
74 are such that, were their values known, our binding 
energies could easily be off by perhaps a few tens of mil- 
livolts. We expect, however, that the general features 
found here should still pertain, namely a single bound 
state whose binding energy is maximized at several tens 
of millivolts along the K-M edge in the basal Brillouin 
zone. 
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